; / 

NASA Contractor Report 191555 7 

ICASE Report No. 93-78 


ICASE 


Years of 
Excellence 


SHAPE OPTIMIZATION GOVERNED BY THE EULER EQUATIONS 
USING AN ADJOINT METHOD 


Angelo Iollo 

OJ 

o 

o 


CO 

Manuel D. Salas 

0 

PsJ 

1 

CO 

ro 

NO 

r-4 

CO 

Shlomo Ta’asan 

1 

o 

z 

o 

c 

D) 

o 

H 

o 


o 

N 

fO 

O 


NASA Contract No. NAS 1 - 1 9480 
November 1993 

Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, Virginia 23681-0001 

Operated by the Universities Space Research Association 



National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23681-0001 


a 

U U Q 

-J C 
X 
LU ►— 

OJ 
Uj T 
X CL 
H *- 
Z ^ 
>- *-* ro 
uj no O 
a 

< o o ^ 

x LU < UJ 
i/) Z Kn 
CC Z < 
14 J < o 

vo a o ^ 

vO O Z 

vO *— * 

H 7 Oi JJ 

c> a z> u 

r-< *"• O 
I I— on CL 
ct < Z <D 
O ^ CD acT 

I k-i ♦— t 

< s: i- — 

CO H <J (t 

< 1— D C 

zaa- 

v □ iii U- 





SHAPE OPTIMIZATION GOVERNED BY THE EULER EQUATIONS 

USING AN ADJOINT METHOD 


Angelo Iollo 1 

Dipartimento di Ingegneria Aeronautica e Spaziale, Politecnico di Torino 
Institute for Computer Applications in Science and Engineering 


Manuel D. Salas 
NASA Langley Research Center 


Shlomo Ta’asan 2 

The Weizman Institute of Science 
Institute for Computer Applications in Science and Engineering 


ABSTRACT 

In this paper we discuss a numerical approach for the treatment of optimal shape 
problems governed by the Euler equations. In particular, we focus on flows with embedded 
shocks. We consider a very simple problem: the design of a quasi-one-dimensional Laval 
nozzle. We introduce a cost function and a set of Lagrange multipliers to achieve the 
minimum. The nature of the resulting costate equations is discussed. A theoretical difficulty 
that arises for cases with embedded shocks is pointed out and solved. Finally, some results 
are given to illustrate the effectiveness of the method. 
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1 Introduction 


The physical structure of the complex flows that occur in aerodynamic design can be 
predicted by reliable numerical simulations. On the other hand, the increasing capability 
of computers to perform even larger calculations radically changes the aerodynamic design 
process. Indeed, for engineering purposes, if one can predict performance, it is fundamental 
to know which modification of an aerodynamic configuration improves performance. 

This question has, of course, been addressed long before the advent of computers which 
has led to a broad category of methods known as inverse design. An exhaustive account 
of the historical development of these approaches is given in [1]. Here it is suffice to say 
that these methods, pioneered by Lighthill [6], require knowledge of a desirable pressure 
or velocity distribution. The adequacy of the distribution chosen is dependent on the 
experience of the designer; the resulting shape strongly depends on this choice. An original 
example of such an approach is found in [3] 

The numerical approach that we will use in this paper, lift the dependence on heuristic 
choices of the desirable distribution, allowing the imposition of constrains to be satisfied 
by the solution found. The numerical simulation of the flow and a numerical optimization 
code are coupled. The optimization code calculates the best perturbation to the geometry 
to decrease a cost function. The geometry itself is described by a set of shape coefficients. 

The optimization code can be devised in one of several ways. A common approach is 
to perturb one shape coefficient at a time and compute the derivative of the cost function 
with respect to this coefficient as a finite difference. Although such codes are simple to 
devise, the procedure is costly and can introduce large errors. In a further evolution of this 
approach, an equation is first calculated for the derivative of the cost function with respect 
to the shape coefficient and then solved numerically. An equation must be solved for each 
shape coefficient. A recent application of this method to a two-dimensional supersonic 
problem is found in [2]. 

The approach presented in this paper is a classical optimal control method. We will 
introduce costate variables (Lagrange multiplier) to achieve a minimum. This method has 
been successfully applied in the design of an airfoil in a subsonic potential flow [10]. 

Here, we consider a flow with embedded shocks where the governing equations are the 
Euler equations. We show how to derive an analytical expression of the cost function 
derivatives with respect to the shape coefficients. For this purpose, we solve only one set 
of costate equations. In [9], some difficulties are outlined, and a method to avoid them 
is proposed. In the present approach an optimal shape can be found for problems with 
embedded shocks, without additional complications. A careful examination of the structure 
of the costate equations suggests a method for integrating them with a robust algorithm 
developed for fluid dynamics purposes. 

A comparison of optimization-based approaches for aerodynamics design problems is 
given in [7], although some results for flows with embedded shocks have been questioned 
recently in [9]. 


1 



2 Problem Statement 


We investigate a new method for aerodynamic design and optimization based on the Eu- 
ler equations. In order to demonstrate the ideas of the method, a very simple problem 
is considered: the design of a Laval nozzle assuming inviscid, quasi-one-dimensional flow. 
The optimization problem consists of finding a set of design variables (in this case shape 
parameters) that minimize some cost function, e.g., a desired pressure or velocity distribu- 
tion along the centerline and possibly requiring that some side constraints be satisfied. We 
attack the optimization problem with the adjoint method. The adjoint method introduces 
a new set of equations and unknowns that are solved together with the flow-field equations. 
To better understand how to obtain a solution of the adjoint equations, the properties of 
these equations are discussed. 

Let the extension of the Laval nozzle in physical space be fi = [0,/]. An energylike 
functional denoted by £, is the cost function we want to satisfy. An optimal shape of the 
nozzle is reached when we meet the necessary conditions for a minimum of £. Let 

£ = {p-p*fdx 

where p* is a target pressure distribution along the abscissa x and p is the pressure field 
for the present geometry. The choice of the functional does not affect the generality of the 
method once the dependence of £ with respect to the flow-field variables is determined. 
Nevertheless, the choice of the functional itself (e.g., \p — p*\ instead of (p — p*) 2 ) can affect 
the performance of the optimization algorithm by changing the curvature of the energylike 
surface. 

In the case of a quasi-one-dimensional flow, if we denote 


p = density 
u = velocity 
e = specific total energy 
p = pressure 
a = speed of sound 
h = height of the channel 
7 = specific heat ratio 


then the Euler equations for unsteady flows reduce to 


where 


U t + F x + Q = 0 


U = 



F = 



Q = 


pup \ 

pu 2 P 

u(pe -(- p)p ) 


( 1 ) 


2 



/? = h x /h, and p = np( 2e — u 2 ). In the following derivation, we use the homogeneity 


property 

F = A(U)U 

( 2 ) 

where 

dF 

( 3 ) 


- A < u) 

and 




/ 0 1 0 \ 
A(U) = I |( 7 ~3)u 2 (3 - 7 )u 2k 1 

\ 2 ku 3 — 7 ue/p ~fe/p — 3 ku 2 7 u } 


The source term Q can be written to display its dependence on U; in fact, the multi- 
plication can be carried out to show that 


Q = S(U)U 

where 


( 4 ) 




0 


0 

— u 2 

\ 2 ku 3 — 'yue/p 


1 

2u 

7 e/p — 3 ku 2 



( 5 ) 


The first and the third row of A(U) are proportional to the first and third row of S(U), 
respectively. Furthermore because d(pu 2 ) — — u 2 • dp-\-2u- d(pu) and pu 2 = — u 2 -p + 2u-pu, 
it follows that S(U) = 5Q/5U. 

We refer to the solution of the above equations as the analysis problem. The bound- 
ary conditions for this problem must be chosen such that the problem is well posed. In 
particular, we will consider the inlet flow with a constant total pressure and entropy, i.e., 
p° tn = p in ( 1 + kM = Constant, Pi 7l /p] n — Constant. At the outlet, if the flow is 
subsonic, the static pressure is fixed as p oui = Constant. 


3 Adjoint Formulation for a Shockless Case 

The design problem can be thought of as a search for a minimum of a functional under 
constrains. Let 

£(«,) = £+ I A t (F x + Q )dx (6) 

Jn 

where A r is an arbitrary vector with components (Ai(x), A 2 (x), A 3 (x)), f! is the domain 
[0, /], and a, are shape coefficients that define the geometry of the nozzle, for example by 
k(a t ,x) = J2i a ifi( x ) fi( x ) a generic function of x. Since in the steady state the 

Euler equations must be satisfied everywhere in the domain, the functionals £ and C are 
identical. If we increase the shape coefficients by eai, then the latter functional will change 
by an amount, say, e8C. The other quantities will_change in the same way: (i =7 ft + £/3, 
p => p + ep, and U(x) =7 U(z) + £U(x), where U = (p, pu, pe) T . By using eqs.(3) and 
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(5), we obtain also F(x) =$- F(x) + eA(U)U(x) and Q(x) => Q(x) + eS(U)U(x). If we 
substitute these relations into eq.(6) and retain only the first-order terms, we obtain 


8£= f l p(p- p *)dx+ /V{[A(U)U] r + S(U)U }dx+ /3A S f U ) U ( ix 
Jo Jo Jo p 


(7) 


In the above equation, /? = ( h x h — hh x )/h 2 . With the notation c,(x) = {hfi x — h x fi)/h 2 , 
the last integral of eq.(7) can be written as 


Wo'- 


c,A t S(U)U 


0 


dx 


( 8 ) 


in which the substitution (3 = <5,c, is made. 

Let us integrate by parts the second integral in eq.(7), with A = A(U) and S = S(U) 
we obtain 


/ ‘ A T [(AU) I + SU]dx = [A t AU]|, - AlAUdx -f [‘ A T SU dx 
Jo Jo Jo 


(9) 


The first term on the right-hand side of the above equation drops for a suitable choice of 
A at the boundaries. The boundary conditions for U are complementary to those for A, 
in the sense that A T AU = 0 yields a homogeneous linear system in A whose rank depends 
on the number of boundary conditions for U. 

For this test case, at the inlet we have 

p° n = Constant 

which implies that 

p = —puu 

If the entropy is fixed, then the specific total enthalpy is also constant; therefore, 

1 


7P 


We conclude that 


(t -Op 2 


[AU]o = [pu,upu, (- 


+ -tr = Constant 


IP 


(7-1 )P 2 

Hence, the suitable choice at the inlet for A is 


+ \ u2 )p u f 


Aj -f- 11X2 + 


IP 


1 2 

+ w 


A3 — 0 


.(7-1 )P 2 

At the outlet for a subsonic case with a given p ou u we obtain 


( 10 ) 


[AU]; = {pu,2upu — u 2 p, 


IP 


3 2 

+ w 


(7 — l)p 2 



7 P , 2 

pu — u 

[h-l)l. +u J 


p} : 
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which leads to 


Aj + 2,11X2 + 

u\ 2 + 


IP 


+ ^ 


L(t-i)p 2 


ip 


+ u' 


A 3 — 0 

a 3 = 0 


( 11 ) 


L(7- l )p 

If the outflow is supersonic, then no boundary conditions are required for U; therefore, 

A = 0 (12) 


identically at the outlet. 

If we take boundary conditions (10) and (11) and eq.(9) into account, eq.(7) can be 
rewritten as 

SC = jf U t [-A t A i + S t A + ^j(p-p*)]dx + X) 5 * J o C ‘ A / g SU (jx 0 3 ) 


where dp/d U = — 2u, 2). If we select A such that 


A t A, + S t A + -p") =0 


(14) 


eq.(13) becomes 


c r _sr djC ~ 

6C ~V da*' 




^A r SU 

0 


dx 


in which we recognize 


OC c t A r SU f 

— — = / dx 

oai Jo p 


(15) 


Suppose that the flow field is known, such that all the variables, dependent upon U, are 
fixed. If we solve eq.(14) with the appropriate boundary conditions for A and substitute 
the results in eq.(15), then we obtain a formulation for the gradient of the Lagrangian. The 
problem then is reduced to finding the solution of a linear system of equations in A with 
homogeneous boundary conditions. Note that A = 0 is a solution if p = p*, which is a 
sufficient condition for SC — 0. If at the minimum p ^ p* although the integral in eq.(15) is 
0, then in general A / 0. The discussion thus far leads to some basic questions about the 
well posedness of the eq.(14) with the boundary conditions given by eqs.(10), (11) or (12) 
and about the existence and uniqueness of the solution. To actually determine a solution 
of this system, examine eq.(14), embedded in time as 


± A, 


A r A x + S r A + q^j(p 


P ’)= o 


( 16 ) 


in which we must choose for the proper sign for the time derivative. 



The matrix A T has the familiar eigenvalues u — a, u and u + a. At the boundaries, if we 
choose the positive sign in eq.(16), for each boundary condition there will be an incoming 
characteristic, such that the problem is well posed. Another motivation for this choice 
is that if we add to eq.(l) the diffusive term, when integrated by parts twice, it would 
append to eq.(14) a second derivative term in x with the same sign this term had in the 
flow equations. Therefore, if the negative sign in eq.(16) is adopted, then an ill-posed heat 
equation would result. In conclusion, note that the costate equation 

At - A r A x + S r A + | *j(p — p* ) = 0 (17) 

has an upside-down characteristic pattern with respect to the time-dependent flow equation. 
If we consider a transonic nozzle, in the throat area, the eigenvalue u — a goes continuously 
through zero; the flow undergoes an expansion through a transonic fan. On the other hand, 
the behavior of the same family of characteristics for eq.(17) shows a shocklike pattern. The 
two characteristic patterns are illustrated in figs. 1(a) and 1(b). 

For a reason that will be clear later, some ambiguous jump conditions for this shock 
will be derived. First consider a simple equation of the kind <!>* — x$ x = 0, with $ = 4>(x), 
x G [—1,1], and boundary conditions 4>(— 1) = $/, 4>(1) = $ r . The characteristics at the 
boundaries show that this problem is well posed and independent of the initial conditions, 
for a large <, the solution will be a step function <I>(x) = for x € [— 1,0[ and $(x) = <l> r 
for x G]0, 1], Note that the jump at zero is solely determined by the boundary conditions 
and that the steady solution will be reached for t — ► oo because of the nearly vertical 
characteristics next to x = 0. 

The structure of the solution to this equation is similar to that underlying eq.(17) for 
which the analysis hides somehow this behavior. A solution to eq.(14) can be written in 
the form 

A(x) = C(x) + [A ]H(x - xth) 

where C(x) G C 1 , [A] is the jump at the shock, H(x — x t h) is the Heaviside function, and x t h 
is the location of the shock for A (i.e., the throat of the nozzle). If we substitute in eq.(14), 
because £( 2 ) (which is the Dirac measure of x) is the derivative of H(x) with respect to x , 
we obtain at x = x t h 

A r [A]£(x - x t h) = 0 

in which all the negligible terms have been dropped. Note that the presence of source terms 
in eq.(14) does not affect the derivation of the jump conditions. 

Finally, because at x = x t k det A = 0, we obtain a nontrivial solution for the system 

A t [A] = 0 (18) 

which yields only two jump conditions; the third jump condition depends on the boundary 
conditions. In this sense this problem has ambiguous jump conditions. 

If the only solution of the homogeneous problem associated with eq.(14), i.e. 

-A t A x + S t A = 0 (19) 
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and its boundary conditions, is A = 0, then thk solution with a nonhomogeneous source 
term is unique, since eq.(14) is linear. Let us consider a subsonic case, with det A / 0 
everywhere in the domain. A general solution of eq.(19) can be written as 

A = A 

and, together with boundary conditions (10) and (11), this implies A = 0 on the domain. 
In the transonic case, since detA = 0 at x = x t h , we split the problem into two domains, 
such that in the subsets [0, x t h[ and ]x t h,l], detA ^ 0. The solution will be 

A su t, = A ' 0 j; ,h ^ T r‘s^ for x e 

and 

ff ft (A T ) _ 1 S T dx r -■] ii 

A 5uper = A 0 e Jx ‘* for x e]x tk , /]. 

Again, if we account for the inlet boundary conditions (10), the jump conditions from 
eq.(18), and the outlet boundary conditions (12), then the solution is A = 0. 

In summary, we have derived an analytic formulation for the gradient of the Lagrangian 
in eq.(6) with respect to the geometry. Furthermore, we have shown that this representation 
is unique in the sense discussed above. 


4 Costate Equations for a Shock Case 


Until now, we have limited our investigation to shockless nozzles to avoid certain difficulties 
that we will discuss here. One problem is that eq.(l) and, therefore, eq.(16) are not defined 
at the shock. This problem is overcome by extending the solution space of U(x) to a set of 
generalized functions, such that eq.(l) will reduce to the Rankine-Hugoniot jumps at the 
shock. A more subtle shortcoming is better understood with the aid the following example. 
Consider a simple equation of the kind + (//(#) — 1 /2)$ x = 0 that is defined, for example, 
on fl = [—1,1]. The characteristics pattern (fig. 2), shows the necessity of some boundary 
conditions on both sides of the the discontinuity to ensure the existence of a steady-state 
solution, regardless of the boundary conditions at the ends of the domain f?. 

Now, eq.(17) can be rewritten to reflect its characteristics pattern 


where 


dp 


P t A, - DP T A r + P r S T A + P 7 -^-(p - p*) = 0 

a U 


P = 


/ 1 1 1 

u — a u u + a 

\(e + p)-ua \u 2 ( e + p) + ua 


( 20 ) 


is the matrix of the right eigenvectors of A(U), and 
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/ u — a 0 0 \ 

D = 0 u 0 

\ 0 0 u + a J 


At the shock wave, the characteristic that corresponds to the eigenvalue u — a undergoes 
a jump in speed of the kind described in the example. If this point is considered inside 
the domain of calculation, then some condition is needed to update the solution in time. 
A boundary condition is needed at this point to continue the calculation to the left of the 
shock wave. We cannot try to derive some boundary conditions for this point, as explained 
for the inlet and the outlet. No speculation about the perturbation U at the shock is 
possible because a perturbation that is solely dependent on the shape coefficients a t and 
the flow equations would be chosen arbitrarily. No other constrains exist. For example, 
the application of the Rankine-Hugoniot jumps to U on each side of the shock would be 
equivalent to assuming that the shock does not change position regardless of the value of 
Oi. The integrals in eq.(7) are split in two, and the integration is carried between [0, x s h[ 
and ]x s h, l ] as 

SjC = SjC\ T &C,2i 

where 


sc, = [a t au]5"* + J*’ h U T [— A t A x + S t A + — P*)]^ 

t*.h C, A T SU 
z^ a < / 

i Jo 


and 


+ 




~dx 


( 21 ) 


SC 2 = 


[A T AU]i.„ + /' U T (- A r A, + S t A + &(p - p'))dx 

+ / 

j X 


X 9 h 

1 CiA r SU 




( 22 ) 


A suitable choice for the Lagrange multiplier is to take A = 0 on both sides of the shock 
such that A is continuous. This selection frees us from imposing a condition on U, because 
the addendum [A T AU] in eqs.(21) and (22) drops anyway at the shock. At the shock, we 
will have three characteristics that deliver the information A = 0 to the left domain, and 
one that delivers it to the right. 

We have examined also another possible interpretation of eqs.(21) and (22). Assume 
that for some perturbation of the shape coefficients the shock properties do not change, 
i.e., the jump in total pressure across it, is essentially constant, which will be the case 
for a weak shock. If we assume that the total enthalpy is constant, then the field to 
parformulationthe right of the shock can be regarded as a subsonic nozzle governed by 
eq.(17) with boundary conditions (10) and (11). To the left of the shock, the flow behaves 
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as a supersonic outlet, where we must impose the proper boundary condition, eq.(12), as 
discussed earlier. With this approach, 8 C\ — 0 and 8 C 2 = 0 independently at the minimum. 

The two approaches presented to handle the shock are not significantly different. In 
the numerical calculations that follow, we have implemented the second set of boundary 
conditions. 


5 Numerical Approach 

The flow-field solution is obtained by introducing a discrete grid (x n , i*) = (x 0 + nAx, t 0 + 
A tk)i where Ax is constant and A tk changes to satisfy the CFL condition. The con- 
servative variables U(x) are computed at the cell centers and integrated in time with a 
three-stage Runge-Kutta scheme, as explained in [5]. In this implementation, we interpo- 
late U(x) to the cell faces by using characteristic differences and a minmod limiter. The 
flux derivative in eq.(l) is then computed using an approximate Riemann solver. See [4]. 
Away from discontinuities, the scheme is second order accurate. 

Depending on the case considered, the solutions of the costate eq.(14), are sought as 
steady results of eq.(17), with boundary conditions applied as explained in the previous 
section. Although eq.(17) is linear, it presents some numerical difficulties because of the 
characteristics pattern at the throat and at the shock location. In particular, consider 
eq.(20) discretized over the same uniform grid of the analysis problem with spacing Ax. If 
we denote by A(-) v the finite increments of the function (•) with respect to the superscripted 
variable, we have, at each grid point 

p T ^_ Dp T ^ + p TsTA + p T ^ (p _ p , ) = 0 (23) 

Define the local increment AW = P r AA; with this notation eq.(23) becomes 


AW 1 

At 


D 


AW 1 T T T dp , 

__ + pr s r A + p ,_> 


P') = 0. 


(24) 


This equation describes the signals that propagate along the characteristics; therefore, 
the increment AW 1 , is one-sided depending on the sign of the corresponding propagation 
speed. Note that for this equation it would be impossible to use a conservative scheme 
since no conservation law exists to satisfy. The integration in time is made by explicit time 
stepping. The scheme is first-order accurate. 

At the throat of the nozzle, u — a = 0 and ft = 0; hence, the first row of eq.(24) reduces 
to 


AW' 

At 


+ pT dv {p - p '> = 0 


The eigenvalue u — a has been shown to vanish at the throat. For a grid point in the 
neighborhood of the throat, this singularity can lead to unbounded grows for Lambda, 
depending on the nature of the source term in the above equation. Because the other two 
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characteristics are nonzero at the throat, this error can degrade the entire calculation. To 
avoid this problem, we approximate u — a with its value at the neighboring point on the 
side from which the vanishing characteristic propagates. 

Since with a Godunov-type solver the shock is resolved with three grid points, we must 
decide on which of these to impose the boundary conditions for A. We must consider that 
the middle point of the three cells on which the shock is solved, is almost sonic; if the 
boundary conditions for A at a subsonic inlet (eq.(10)) were imposed at this grid point, 
then the convergence rate to the steady solution would be considerably slower. For this 
reason we impose the condition for supersonic outlet A = 0, on the middle grid point, 
eliminating it from of the computation. 

Another remark should be made in regard to the order of magnitude of the residu- 
als of eq.(23), for which we can consider the solution steady. Close to the minimum, the 
gradients in eq.(15) are almost zero; nevertheless the optimization algorithm requires a 
careful computation of these values, such that in order to consider the time-dependent so- 
lution converged, the residuals must be some orders of magnitude smaller than the gradient 
components. 

In the results that follow, we used a representation of the nozzle geometry defined by 
h = a^X + a-ijX + 03, where X = x + 10 -3 ; this representation allows two independent 
design variables because — h x /h. 

In this work, we do not address the methods to accelerate the numerical scheme to 
obtain an optimal shape; the strategy used to achieve the minimum of the functional is 
straightforward: 

1 . Start with a first guess for the shape coefficients. 

2. Solve the flow equation. 

3. Solve the costate equation with the values computed in step 2 for the flow field. 

4. Update the shape coefficient with a gradient-based criterion. 

5. Restart the procedure from step 2 until the gradient is zero. 

To update the shape coefficients a BFGS algorithm was used. See [8]. In some cases, 
as we will discuss, we used an inefficient, but robust, algorithm that simply makes a line 
search for a zero of the gradient. 


6 Discussion of the Results and Conclusions 

The values for the functional £, computed with an analytical solution of eq.(l), are shown 
in fig.(3). The numerical values of the analytical solution are computed on the same grid 
presented above; then, the functional is computed by a trapezoidal approximation. The 
discrete functional, which is a result of the trapezoidal integration, shows some disconti- 
nuities and a local minimum that disappears as the number of grid points increases. As 
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the mesh is refined, the number of the discontinuities increases, while the jumps become 
smaller. 

Note that if the dependence of the geometry on the shape coefficients is smooth, then 
the functional 



is always defined with the assumption that dpji 9cv, is defined everywhere except at a finite 
number of points. This assumption is reasonable because the solution of eq.(l) changes 
smoothly with the geometry. For this reason, even if no certainty exists that the solution 
depends monotonically on the shape coefficients, this behavior can be the interpreted in 
this way. Suppose that the integrand of the functional can be represented by a simple 
rectangular function. If in attaining the minimum the area under the curve decreases and its 
“height” increases, the functional will eventually increase before the edges of the rectangle 
pass another grid point, because the mesh resolution is not sufficient. The functional will 
exhibit a local minimum and a subsequent discontinuity. 

A method that derives a formulation for the gradient of £ from a discrete approximation 
of the functional will obtain meaningless solutions as a result of the discontinuities of the 
discrete functional, such that no optimization algorithm alone could anyway get to the 
minimum. In the present formulation, an approximate representation of the analytical 
gradient of the functional is derived. For this reason, the approximation of the analytical 
gradient will be, at most, affected by discontinuities due to the discretization and will be 
always monotonic (if the analytical functional does not change curvature) and bounded. See 
fig. 4. In figs. 5 and 6, we present two sets of results in which the target pressure distribution 
was generated with the same h(x) that was used in the optimization procedure. In fig. 5, 
the target pressure distribution is obtained starting from a subsonic first guess. This result 
shows the effectiveness of the method. Fig. 6 shows that we can achieve the optimum from 
both sides of the discontinuity. 

In general, when p*(x) is fixed, the minimum of £ is reached for different values of the 
shape coefficient a,. These different values depend on whether one considers the analytical 
or the discretized functional, even if the results are converged on the grid. If the target 
solution can be attained exactly, then both values coincide. The gradient calculated after 
the proposed derivation, will still depend on the discretization through the nonhomogeneous 
term in eq.(14). In fig. 7, we show, for a case in which p* cannot be reached exactly, that 
the distance between the two minima becomes approximately half when the grid resolution 
is doubled. This result supports the hypothesis that the minimum calculated through the 
analytical gradient will indefinitely approach to the actual minimum as the grid is refined. 

In conclusion, a method has been presented to calculate the gradient components of a 
generic functional, in which (regardless of the number of the shape coefficients) only one 
linear costate equation must be solved. The minimum computed in this way differs from 
the minimum of the discrete functional; however these minima indefinitely approach as the 
grid is refined. 
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(a) 


(b) 




Figure 1: Characteristic patterns, (a) Transonic expansion, (b) Correspondent shocklike 
structure for costate equations. 



Figure 2: Characteristic pattern for + [ H(x ) — 1/2]4> X = 0. Boundary conditions are 
needed on both sides of discontinuity. 
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(a) (b) 




Figure 3: Analytical solution calculated with shape function h(x) = a(x 3 — x 2 ) + 1.05; 
functional has been calculated with trapezoidal approximation, (a) Solution distributed 
over 20 grid points, (b) Solution distributed over 80 grid points. 
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(a) (b) 



Figure 4: Result shown is obtained with h(x) = a\X + 0.3 j X + 10. For each value of a \ 
flow equation is solved by a Godunov like scheme, and gradient is calculated as proposed 
in this paper. Each time the shock goes through grid point, discretized functional does 
not have a monotonic derivative; gradient has discontinuities, but is still monotonic, (a) 
Functional, (b) Gradient. 
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Pressure 



Figure 5: Target solution is h(x) = 2.bX + 0.3/ X + 5. First guess is h(x) = 2X + 0.52/X + 5. 
Shape function for which minimum is sought is h(x) — ot\X + ci^jX + 5, and optimization 
algorithm used is BFGS. Starting value of the functional is of order 10~ 3 . At minimum it 
is of order 10 -9 . Gradient components are of order 10 -12 at minimum. 
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Pressure 



Figure 6: Target solution is h(x) = 2 . bX + Q.'S/ X-\- 10. First guess is h(x) = 5 A" H- 0.3/ A" + 1 0. 
Shape function for which minimum is sought is h(x) = a\X + a 2 / X + 10, and optimization 
algorithm used is BFC1S. Starting value of functional is of order of 10 -3 . At minimum it is 
of order 10 -9 . Gradient components are of order 10 -12 at minimum. 
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Figure 7: Target solution is h(x) = 2.5X + OA/X + 10, and the shape function for which 
minimum is sought is h[x) = a\X + 0.3/X +5. Divide and conquer procedure was used to 
find zero of gradient, (a) Convergence of gradient to zero, (b) Functional value, at zero of 
gradient, is reduced to half by doubling grid points, (c) Design variable. 
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